Dynamical scaling in branching models for seismicity 
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We propose a branching process based on a dynamical scaling hypothesis relating time and mass. 
In the context of earthquake occurrence, we show that experimental power laws in size and time 
distribution naturally originate solely from this scaling hypothesis. We present a numerical protocol 
able to generate a synthetic catalog with an arbitrary large number of events. The numerical data 
reproduce the hierarchical organization in time and magnitude of experimental inter-event time 
distribution. 

PACS numbers: 02.50.Ey,64.60.Ht,89.75.Da,91.30.Dk 



Stochastic branching models have been used since a 
long time in the description of a large variety of social and 
physical phenomena ranging from biological evolution or 
genealogy [J] , to nuclear or chemical reactions [1] . In the 
last years branching processes have been widely studied 
in seismicity and actually they provide one of the most 
efficient tools for earthquake forecasting 0]. Within this 
approach, one treats seismicity as a marked point pro- 
cess in time {Mi(ti)} where ti is the occurrence time of 
an event with mass (magnitude) Mj and one assumes 
that each event can trigger future ones according to a 
two point conditional rate p{M(t)\M i (t i )). Given a his- 
tory of past events {Mj(ij)}, then, the rate of events of 
magnitude M at time t is given by 



p{M{t)\{Mi(M)}) 



p(M(t)|M i (* i ))+AiP(M) (1) 

i:U<t 

where p is a constant rate of independent sources and 
P{M) their magnitude distribution. In the epidemic type 
aftershock sequences (ETAS) model [4J, widely used in 
seismology, for any couple of events i and j, magnitude 
of event i is independent of previous events 



p{M i {U)\M j {t j )) = P(Mi)g{ti - tfM,) 



(2) 



and for the magnitude distribution P(Mj) and the prop- 
agator g(ti — tj]Mj) one uses three well known experi- 
mental observations: 

i) The magnitude distribution follows an exponential 
law P(M) ~ 10 -bM usually referred as the Gutenberg- 
Richter (GR) law @, where b ~ 1; 

ii) The Omori law Q states that the number of af- 
tershocks n(t) decays in time as n(t) ~ (t + c)~ p with 
p~l; 

iii) The aftershock number is exponentially related to 
the mainshock magnitude. This last law combined with 
the Omori law gives </;/, /,...!/,.: oc 10 aM i (U-tj + c)- p , 
where a ~ b 0]. 

The main statistical properties of the ETAS model 
have been reviewed in a series of papers (see for in- 
stance ref. [§]). Ultraviolet and infrared cut-offs have 



to be necessarily introduced to make the theory con- 
vergent. Whereas a large magnitude cut-off can be ex- 
pected on physical grounds, more questionable is the ex- 
istence of a minimum magnitude Mj n / Q . This problem 
has been removed by a self-similar version of the ETAS 
model recently introduced by Vere- Jones [loj . In this 
model the decoupling (|2|) between magnitude and time 
is still assumed but a multiplicative factor is considered 
p{M i {t i )\M j {t j )) = P(M i )g(t i -t j ;M j )S(M i -M j ) where 
S(Mi - Mj) = W- d \ M *- M i\ introduces magnitude corre- 
lations: for d > large daughter earthquakes tend to oc- 
cur after large mother earthquakes. Saichev and Sornette 
have suggested a physical explanation for this magnitude 
correlation based on faults branching [ll| . 

In this paper we show that the above mentioned magni- 
tude correlations do not need to be introduced via an ad 
hoc term, as S(M, — Mj), but naturally originate from a 
more general scaling relation. More precisely we assume 
that the magnitude difference Mj — Mj fixes a character- 
istic time scale 



km b{M,- Mi ) 



(3) 



so that the conditional rate is magnitude independent 
when time is rescaled by Tjj and k is a constant measured 
in seconds 



p{Mi{U)\Mj{tj)) = F 



U - L 



(4) 



With the only constraint that F(x) can be normalized, 
we recover the statistical features of earthquake occur- 
rence: Omori law, GR law and scaling behaviour of the 
interevent time distribution 12, Tj^. Furthermore, using 
Eq.s (|l|4j) in a numerical code we are able to generate, in 
few hours of CPU time, a synthetic catalog with the same 
number of events as 30 years California catalog. Exper- 
imental and numerical catalogs are found to exhibit the 
same time and magnitude organization. 

In order to verify the existence of magnitude corre- 
lations we consider the ANSS California catalog 14 1 
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FIG. 1: (Color online) Function C(n) versus n for the ANSS 
catalog (o), the restricted catalog (□) and the simulated cat- 
alog (A). C(n) evaluated for five different realizations of 
reshuffled ANSS catalog (dots). 



containing N — 9586 M > 3 earthquakes. We di- 
vide the catalog in Nl = N/L subsets each contain- 
ing L = 125 events and define the quantity 6Mj — 

- (V^)Eti^ representing the 
deviation of the average magnitude in the j-th subset 
with respect to the average over the entire catalog. We 
then calculate the quantity C(n) defined as 



C(n) 



n+l 



I N L - 

j=n+l 



N L -n m 

E 

i=l 



(5) 



where n max — 32 is the maximum "distance" between 
subsets considered. The advantage of correlating aver- 
age quantities SMj is to reduce fluctuations and the sum 
over j further smoothens statistical noise. In absence 
of magnitude correlations, 8Mi and SMi + j are statisti- 
cally independent quantities and C(n) does not depend 
on n and fluctuates around C(n) = 0. This situation 
can be realized by evaluating C(n) after reshuffling the 
magnitudes. Using 100000 realizations of the reshuffled 
catalog, we find that the distribution of C(n) exhibits 
gaussian behaviour centered in zero with standard devi- 
ation a = 0.0001. Figure 1, conversely, shows that C(n) 
computed for the real catalog has a regular trend with 
amplitude several times larger than cr, clearly indicating 
the existence of magnitude correlations. Therefore the 
behaviour of C(n) for the real catalog is a signature of a 
well defined earthquake magnitude organization. 

In order to check that the results of Fig. 1 are not a 
spurious effect of short term aftershock incompleteness 
15| . we use the method proposed by Helmstetter ct al 
7] stating that, after a main shock of magnitude Mm 
at time the completeness level M > 3 is recovered 
only after a time t c = t M + lO^-" 7 - 5 )/ - 75 . We then 
construct a restricted catalog by neglecting all events oc- 
curring within a time interval [t m , tc] after each event 
with M > 1, obtaining a catalog complete for M > 3 



and containing N = 8502 events. The evaluation of C(n) 
(Fig.l) now shows that correlations do not disappear but 
C(n) exhibits again a regular trend with values signifi- 
cantly different than zero. The comparison of C(n) in 
the original and restricted catalog indicates that magni- 
tude correlations are not therefore an artifact of catalog 
incompleteness but should be attributed to a physical ef- 
fect due to earthquake interactions. These interactions 
are introduced in our approach by means of the scaling 
assumption ^ . In this way we are able to reproduce, at 
least at a qualitative level, the experimental behaviour of 
C(n) as can be seen in Fig.l. 

Before discussing the details of our numerical proce- 
dure, we explore the consequences of our scaling as- 
sumption ([3]). We first notice that the total number of 
daughter earthquakes, conditioned to the occurrence of a 
mother earthquake of magnitude Mq at time to, is given 
by dtp(M(t)\M (t )) = kW-^ M - M °^ J °° dxF(x). 
Since F(x) is normalizable, we recover the GR be- 
haviour independently of the specific form of F(x). On 
the basis of this observation, a mother earthquake at 
time to will be distributed according to GR law and 
therefore the occurrence rate of magnitude M trig- 
gered events at time t is then given by p(M,t — t ) = 
Im'T p(M(t)\M (t ))P(Mo)dM which leads to 

lQ -bM io^-^/X*-^) 
p(M,t-t ) oc -r F(z)dz (6) 

[t - to) J W -HM Bup -M)(t-t ) 

We first observe that ultraviolet and infrared cut-offs 
are not necessary anymore. Indeed, in Eq.([S]) one can 
— oo and M sup — > oo obtaining 

Then one recovers, as a direct 



arbitrary set Mi n f — 
p(M,t-t ) oc 
consequence of the only assumption (jU, beside the GR 
also the Omori law that conversely are assumed a pri- 
ori in the ETAS model. The above result, that hold for 
quite arbitrary F(x), suggests that these two fundamen- 
tal laws, generally considered as independent laws in seis- 
micity, can be strictly related to a more general scaling 
behaviour. Furthermore Eq.^ shows that a change in 
M sup or M in f only corresponds to a time rescaling and 
therefore the only effect is on the parameter k in ©■ 

We now construct with our approach a synthetic cat- 
alog to be compared with the experimental one. In a 
numerical protocol one assumes at initial time to = 
a single event of arbitrary magnitude chosen in a fixed 
range [M in f, M sup ). Time is then increased by a unit step 
t = to + 1, a trial magnitude is randomly chosen in the 
interval [Mi n f,M sup ] and Eqs. (|ll4[ ) give the probability 
to have an earthquake in the time window (to, to + 1). 
If this probability is larger than a random number be- 
tween and 1, an earthquake takes place, its magnitude 
and occurrence time are stored and used for the evalu- 
ation of probability for future events. Time is then in- 
creased and in this way one constructs a synthetic catalog 



of N e events. The term /i in Eq.© represents an ad- 
ditional source of earthquakes Poissonian distributed in 
time with a magnitude chosen from the GR distribution 
with b — 0.8. Other choices for P(M) do not sensibly 
affect the statistical properties of the simulated catalog. 

Following this protocol, we generate sequences of 15000 
events using a power law form for F(z) — A/(z x + 7) 
with Mi n f = 1 and M sup — 8. We compute magnitude 
distribution P(M) and intertime distribution D(At, Ml) 
where At is the time distance between successive events 
with magnitude greater than a given threshold M l . Ex- 
tended analysis of experimental catalogs have shown 
12l T3| that the intertime distribution is a fundamental 
quantity to characterize the magnitude and time organi- 
zation of earthquakes. In fact, indicating with Pc(M) 
the cumulative magnitude distribution inside the consid- 
ered region, one observes 



D(At,M t ) = P c (M L )f(P c (M L )At) 



(7) 



where / is a universal function, independent on Ml and 
on the geographical region indicating a well defined hier- 
archical organization of earthquake occurrence [161 ] . The 
numerical distributions are compared with the experi- 
mental data from the ANSS Catalog. For different val- 
ues of A, it is always possible to find a set of parameters 
A, 7, 6, /i such that numerical data reproduce, on average, 
earthquake statistical features both in time and in mag- 
nitude. The parameter k is fixed a posteriori in order to 
obtain the collapse between numerical and experimental 
time. We have also performed simulations for different 
values of Mi n t and M sup obtaining similar results and 
confirming that changes in the magnitude range only pro- 
duce time rescaling. 

In Fig. 2 we plot the experimental and numerical 
D(At, Ml) considering two different values of A (A = 1.2 
and 5) and M L (M L = 1.5 and 2.5). Data for different 
values of the parameters follow a universal curve and the 
same collapse is obtained for other values of A > 1. The 
accordance between experimental and numerical curves 
(inset of Fig. 2) indicates that the hypothesis of dynamical 
scaling is able to reproduce two fundamental properties 
of seismic occurrence, namely the GR law and Eq. ([7]), 
independently of the details of F(z). 

The ETAS model corresponds to a particular choice for 
F(z), i.e. 7 = and A = p a 1. We want to stress the 
important difference due to the presence of a non-zero 7. 
First, the constant 7 removes the problematic need of an 
infrared (ultraviolet) cut-off in the ETAS model, whereas 
in our approach Mi n f and M sup are irrelevant variables. 
Second, the constant 7 gives rise to the observed magni- 
tude correlation. Indeed, for a given mainshock of magni- 
tude Mj at time tj, at each time (U > tj) it is possible to 
define a sufficiently large magnitude difference AM such 
that, if M. } - Mi > AM, we have that z x is negligible 
with respect to 7 and therefore F[(ti — tj)/rij] ~ A/j. 
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FIG. 2: (Color online) The intertime distribution with F(z) 



A/(z + 7), with two values of A 



1.2,5 and M L 



1.5, 2.5. Continuous and broken curve are the experimen- 
tal D(At, Ml) with Ml = 1.5 and Ml = 2.5. For A = 1.2 
(A = 5) we set k = 210sec (k = 420sec), A = 1.4 10" 4 sec _1 
( A = 1.9 hnSec" 1 ), fi = 4 10" 7 (ju = 1.510 -6 ), 7 = 1 
(7 = 0.1) and 6 = 1. In the inset the magnitude distribution 
of the experimental (black line) and numerical catalog with 
A = 1.2 (red o) and A = 5 (green □). 



In other words after a large event, the probability of big 
quakes is raised. 

We have also performed more extensive simulations us- 
ing for F{z) an exponential behaviour 



F(z) = 



A 



1+7 



(8) 



Eq.© states that two events of magnitude Mi and Mj 
are correlated over a characteristic time and become 
independent when ti — tj > Tij. As a consequence, only 
a small fraction of previous events can affect the proba- 
bility of future earthquakes so that, after a certain time, 
Earth crust loses memory of previous seismicity. This 
aspect is perhaps more realistic with respect to the idea, 
contained in a power law correlation, that events are all 
correlated with each other and also gives rise to impor- 
tant implications for seismic forecasting. The construc- 
tion of seismic catalogs, indeed, dates back to about 50 
years, and according to Eq.© one can have good es- 
timates of seismic hazard without considering previous 
seismicity. This is no longer true if one assumes a power 
law time decorrelation. We want also to point out that 
a general state-rate formulation [llj gives rise to corre- 
lations between earthquakes that decay exponentially in 
time. We finally observe, that taking into account only 
a fraction of previous events in the evaluation of con- 
ditional probabilities, the numerical procedure consider- 
ably speeds up. In the case of long range temporal cor- 
relations, CPU time grows with the number of events as 
iVg , whereas in the case of an exponential tail the growth 
is linear in N e . For this reason, assuming the functional 
form ijSJl one can simulate very large sequences of events. 
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FIG. 3: (Color online) The intertime distribution as function 
of AtP c (Mi) obtained using Eq.([5} (black circle O) and com- 
pared with the experimental distributions (red □) for three 
different values of Ml (Ml = 1.5,2.5,3.5, from top to bot- 
tom). We set k = 4.9 W 4 sec, A = 6.1 10~ 5 sec _1 , fx = 2 10" 5 , 
7 = 0.1. In inset (a) collapse of the three curves for different 
Ml following the scaling of Eq.([7} and in inset (b) the exper- 
imental (red) and numerical (black) magnitude distribution. 

In particular for a different choice of parameters, one can 
construct synthetic catalogs containing the same number 
of events (N e = 245000 with M > 1.5) of the experimen- 
tal California Catalog. In Fig. [3] we compare numerical 
and experimental distributions D(At, Ml) for three dif- 
ferent values of Ml- For each value of Ml, the numer- 
ical curve reproduces the experimental data and fulfills 
Eq.([7]) (inset (a) in Fig.(j3j)). We have also evaluated the 
behaviour of C(n) using the numerical catalog finding 
qualitative agreement with experimental results (Fig.l). 
Finally, the numerical magnitude distribution P(M) fits 
very well the experimental one (inset (b) in Fig. ([5])). 

After fixing k, we express numerical time unit in sec- 
onds and we observe that the numerical catalog corre- 
sponds to a period of about 9.9 10 9 sec ~ 30 years. Our 
model is therefore able to construct a synthetic cata- 
log covering about 30 years, containing about the same 
number of events and displaying the same statistical or- 
ganization in magnitude and time of occurrence as real 
California Catalog. The high efficiency of the model in 
reproducing past seismicity indicates that the model is 
a good tool for earthquake forecasting. In fact, given a 
seismic history, Eq.([T]) together with Eq.s(Ul[5|) gives the 
rate of occurrence of magnitude M earthquakes at time t 
inside a considered geographic region. We want to point 
out that similar extended analysis has never been per- 
formed in seismicity. It is, indeed, the first time that a 
stochastic model is able to produce a synthetic catalog 
with all the features of the real experimental intertime 
and magnitude distributions. Recently, Saichev and Sor- 
nette [2l| have calculated the intertime distributions with 



an analytical approach to the ETAS models. Under the 
assumption that an earthquake can have at most one first 
generation aftershock they find for the intertime distri- 
bution a behaviour different from Eq.([7]). 

We finally observe that also spatial organization of 
seismic events reveals some kind of scale invariance 
[TgL Eoj ]. This indicates that also spatial distribu- 
tion originates from a critical behaviour of the Earth 
crust suggesting that a dynamical scaling hypothesis as 
in Eq.Q can also work if one appropriately introduces 
spatial dependencies. In this way it would be possible to 
construct seismic hazard maps. 
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